Assignment #2

hotels <- readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/master/data/2020/2020-02-11/hotels.csv')

When you finish the assignment, remove the # from the options chunk at the top, so that messages and warnings aren’t printed. If you are getting errors in your code, add error = TRUE so that the file knits. I would recommend not removing the # until you are completely finished.

Put it on GitHub!

From now on, GitHub should be part of your routine when doing assignments. I recommend making it part of your process anytime you are working in R, but I’ll make you show it’s part of your process for assignments.

Task: When you are finished with the assignment, post a link below to the GitHub repo for the assignment.

Here is my GitHub repo link: https://github.com/ducminhngo871/Machine-Learninng-Lasso.

Machine Learning review and intro to tidymodels

Read through and follow along with the Machine Learning review with an intro to the tidymodels package posted on the Course Materials page.

Tasks:

1. Read about the hotel booking data, hotels, on the Tidy Tuesday page it came from. There is also a link to an article from the original authors. The outcome we will be predicting is called is_canceled.

• Without doing any analysis, what are some variables you think might be predictive and why?

Without doing any analysis, I think previous_cancellations, previous_bookings_not_canceled and deposit_type might be the three variables that might be predictive. The reason for that is becasue previous_cancellations show the number of previous bookings that were cancelled by the customer prior to the current booking. The same thing can be applied for previous_bookings_not_canceled.

Additionally, for deposit_type, if the room has been charged (Non Refund), it will be more likely that the room will not get cancelled a lot of time.

_ What are some problems that might exist with the data? You might think about how it was collected and who did the collecting.

Some problems might exist with the data are:

• As the data was obtained directly from the hotels’ PMS databases’ servers, the hotel’s workers were the one who complete the data. Therefore, it might be biased based on the view of the hotel’s workers.

  • If we construct a model, what type of conclusions will be able to draw from it?

If we create a model, we can find which variables have the highest impact on the decision of cancelling a booking or not.

2. Create some exploratory plots or table summaries of the variables in the dataset. Be sure to also examine missing values or other interesting values. You may want to adjust the fig.width and fig.height in the code chunk options.

library(tidyverse)         # for reading in data, graphing, and cleaning
library(tidymodels)        # for modeling ... tidily
library(usemodels)         # for suggesting step_XXX() functions
library(glmnet)            # for regularized regression, including LASSO
library(naniar)            # for examining missing values (NAs)
library(lubridate)         # for date manipulation
library(moderndive)        # for King County housing data
library(vip)               # for variable importance plots
library(rmarkdown)         # for paged tables
theme_set(theme_minimal()) # my favorite ggplot2 theme :)

First, I will see the first 5 rows of the data to see its variables:

hotels %>% 
  slice(1:5)

In here, I can see the dataset has 32 variables and the one we will focus on in this assignment is is_canceled: Value indicating if the booking was canceled (1) or not (0).

1. Quick exploration

After seeing the dataset for a brief moment, I will take a quick look at distributions of all the variables to check for anything irregular.

Quantitative variables:

hotels %>% 
  select(where(is.numeric)) %>% 
  pivot_longer(cols = everything(),
               names_to = "variable", 
               values_to = "value") %>% 
  ggplot(aes(x = value)) +
  geom_histogram(bins = 30) +
  facet_wrap(vars(variable), 
             scales = "free")

When I see the histogram plots of the variables, I can see that the data was collected for 3 years with the majority of data comes from 2016 and 2017.

Things I have noticed from the graphs:

• Right skewness in lead_time, stays_in_week_nights –> log transform if using linear regression.

• Many 0’s in children, days_in_waiting_list, adr,adults,is_repeated_guest and required_car_parking_spaces –> create indicator variables of having that feature vs. not, ie. a variable called children where a 0 indicates no basement (children = 0) and a 1 indicates a basement (children > 0).

For the categorical variables:In this dataset, there are no categorical variables.

Additionally, we might want to see the number of cancelled booking versus the number of booking not being cancelled.

hotels %>% 
  count(is_canceled)

From here, we can see that the number of bookings that got canceled is close to the number of bookings that did not get canceled. The number of data is sufficient enough for us to do any machine learning models.

3. First, we will do a couple things to get the data ready.

• I did the following for you: made outcome a factor (needs to be that way for logistic regression), made all character variables factoors, removed the year variable and some reservation status variables, and removed cases with missing values (not NULLs but true missing values).

• You need to split the data into a training and test set, stratifying on the outcome variable, is_canceled. Since we have a lot of data, split the data 50/50 between training and test. I have already set.seed() for you. Be sure to use hotels_mod in the splitting.

hotels_mod <- hotels %>% 
  mutate(is_canceled = as.factor(is_canceled)) %>% 
  mutate(across(where(is.character), as.factor)) %>% 
  select(-arrival_date_year,
         -reservation_status,
         -reservation_status_date) %>% 
  add_n_miss() %>% 
  filter(n_miss_all == 0) %>% 
  select(-n_miss_all)

set.seed(494)

# Randomly assigns 75% of the data to training.
hotels_split <- initial_split(hotels_mod, 
                             prop = .5)
hotels_split

## <Analysis/Assess/Total>
## <59693/59693/119386>

#<training/testing/total>

hotels_training <- training(hotels_split)
hotels_testing <- testing(hotels_split)

4. In this next step, we are going to do the pre-processing. Usually, I won’t tell you exactly what to do here, but for your first exercise, I’ll tell you the steps.

• Set up the recipe with is_canceled as the outcome and all other variables as predictors (HINT: ~.).
  
• Use a step_XXX() function or functions (I think there are other ways to do this, but I found step_mutate_at() easiest) to create some indicator variables for the following variables: children, babies, and previous_cancellations. So, the new variable should be a 1 if the original is more than 0 and 0 otherwise. Make sure you do this in a way that accounts for values that may be larger than any we see in the dataset.
  
• For the agent and company variables, make new indicator variables that are 1 if they have a value of NULL and 0 otherwise. I also used step_mutate_at() for this, but there’s more ways you could do it.
• Use fct_lump_n() inside step_mutate() to lump together countries that aren’t in the top 5 most occurring.
• If you used new names for some of the new variables you created, then remove any variables that are no longer needed.
• Use step_normalize() to center and scale all the non-categorical predictor variables. (Do this BEFORE creating dummy variables. When I tried to do it after, I ran into an error - I’m still investigating why.)
• Create dummy variables for all factors/categorical predictor variables (make sure you have -all_outcomes() in this part!!).
  
• Use the prep() and juice() functions to apply the steps to the training data just to check that everything went as planned.

hotel_recipe <- recipe(is_canceled ~ .,
                       data = hotels_training) %>% 
  step_mutate_at(children, babies, previous_cancellations,
                 fn = ~ as.numeric(. > 0)) %>%
  step_mutate_at(agent, company,
                 fn = ~ as.numeric(. == "NULL")) %>%
  step_mutate(country, 
              country_grp = fct_lump_n(country, n = 5)) %>% 
  step_rm(country) %>% 
  step_normalize(all_predictors(),
                 -all_nominal(),
                 -all_outcomes()) %>%
  step_dummy(all_nominal(),
             -all_outcomes())

hotel_recipe %>% 
  prep(hotels_training) %>%
  # using bake(new_data = NULL) gives same result as juice()
  # bake(new_data = NULL)
  juice() 

5. In this step we will set up a LASSO model and workflow.

• In general, why would we want to use LASSO instead of regular logistic regression? (HINT: think about what happens to the coefficients).

• Lasso regularizations are also known as ‘shrinkage’ methods, because they reduce or shrink the coefficients in the resulting regression. This reduces the variance in the model: as input variables are changed, the model’s prediction changes less than it would have without the regularization. While using LASSO model, we can avoid overfit by reducing the variance of a model.

• Meanwhile, logistic regression does not shrink the variables, so we will have a lot of variables than needed.

• Define the model type, set the engine, set the penalty argument to tune() as a placeholder, and set the mode.
  

• Create a workflow with the recipe and model.

hotel_lasso_mod <- 
  # Define a lasso model 
  logistic_reg(mixture = 1) %>% 
  # Set the engine to "glmnet" 
  set_engine("glmnet") %>% 
  # The parameters we will tune.
  set_args(penalty = tune()) %>% 
  # Use "regression"
  set_mode("classification")

## Set up the workflow: 
hotel_lasso_wf <- 
  # Set up the workflow
  workflow() %>% 
  # Add the recipe
  add_recipe(hotel_recipe) %>% 
  # Add the modeling
  add_model(hotel_lasso_mod)

hotel_lasso_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 6 Recipe Steps
## 
## • step_mutate_at()
## • step_mutate_at()
## • step_mutate()
## • step_rm()
## • step_normalize()
## • step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = tune()
##   mixture = 1
## 
## Computational engine: glmnet

6. In this step, we’ll tune the model and fit the model using the best tuning parameter to the entire training dataset.

• Create a 5-fold cross-validation sample. We’ll use this later. I have set the seed for you.
  
• Use the grid_regular() function to create a grid of 10 potential penalty parameters (we’re keeping this sort of small because the dataset is pretty large). Use that with the 5-fold cv data to tune the model.
  
• Use the tune_grid() function to fit the models with different tuning parameters to the different cross-validation sets.
  
• Use the collect_metrics() function to collect all the metrics from the previous step and create a plot with the accuracy on the y-axis and the penalty term on the x-axis. Put the x-axis on the log scale.
  
• Use the select_best() function to find the best tuning parameter, fit the model using that tuning parameter to the entire training set (HINT: finalize_workflow() and fit()), and display the model results using pull_workflow_fit() and tidy(). Are there some variables with coefficients of 0?

set.seed(494) # for reproducibility

#Split the data into 5-fold: 
hotel_cv <- vfold_cv(hotels_training, v = 5)

# potential penalty parameters
penalty_grid <- grid_regular(penalty(),
                             levels = 10)

hotel_lasso_tune <- 
  hotel_lasso_wf %>% 
  tune_grid(
    resamples = hotel_cv,
    grid = penalty_grid
  )

collect_metrics(hotel_lasso_tune, summarize = TRUE)

### How can you select only "accuracy" not roc_curve. 
collect_metrics(hotel_lasso_tune) %>%  
  filter(.metric == "accuracy") %>% 
  ggplot(aes(x = log(penalty), y = mean)) + geom_point()

hotel_lasso_tune %>% 
  show_best(metric = "accuracy")

# Best tuning parameter by accuracy
best_param <- hotel_lasso_tune %>% 
  select_best(metric = "accuracy")
best_param

hotel_lasso_final_wf <- hotel_lasso_wf %>% 
  finalize_workflow(best_param)
hotel_lasso_final_wf

## ══ Workflow ════════════════════════════════════════════════════════════════════
## Preprocessor: Recipe
## Model: logistic_reg()
## 
## ── Preprocessor ────────────────────────────────────────────────────────────────
## 6 Recipe Steps
## 
## • step_mutate_at()
## • step_mutate_at()
## • step_mutate()
## • step_rm()
## • step_normalize()
## • step_dummy()
## 
## ── Model ───────────────────────────────────────────────────────────────────────
## Logistic Regression Model Specification (classification)
## 
## Main Arguments:
##   penalty = 1e-10
##   mixture = 1
## 
## Computational engine: glmnet

Now we could fit this to the training data and look at the resulting model.

hotel_lasso_final_mod <- hotel_lasso_final_wf %>% 
  fit(data = hotels_training)

hotel_lasso_final_mod %>% 
  pull_workflow_fit() %>% 
  tidy() 

In here, I can see a few terms have coefficients of 0 (such as arrival_date_month_September, market_segment_Undefined, distribution_channel_Undefined)

7. Now that we have a model, let’s evaluate it a bit more. All we have looked at so far is the cross-validated accuracy from the previous step.

• Create a variable importance graph. Which variables show up as the most important? Are you surprised?

# Visualize variable importance
hotel_lasso_final_mod %>% 
  pull_workflow_fit() %>% 
  vip()

In here, the reserved_room_type_P comes up as the most important features, followed by assigned_room_type_I and deposit_type_Non.Refund. For me personally, I am pretty surprised about the variable importance. Since I thought that deposit type should be the most important variables, however, I am fairly surprise that it is only a top 3 categories. Moreover, I thought previous_cancellations might be an important variable, however, I did not see that in the model.

• Use the last_fit() function to fit the final model and then apply it to the testing data. Report the metrics from the testing data using the collect_metrics() function. How do they compare to the cross-validated metrics?

# Fit model with best tuning parameter(s) to training data and apply to test data
hotel_lasso_test <- hotel_lasso_final_wf %>% 
  last_fit(hotels_split)

# Metrics for model applied to test data
hotel_lasso_test %>% 
  collect_metrics()

# This one is not totally right. 
collect_metrics(hotel_lasso_tune) %>%
  filter(.metric == 'accuracy')

In here, when we look at the accuracy rate, the testing data’s accuracy is fairly consistent with the training data. For the testing data, it is estimated that the accuracy is 0.8138140, slightly lower than the training data (0.8161091). However, the number is fairly consistent.

• Use the collect_predictions() function to find the predicted probabilities and classes for the test data. Save this to a new dataset called preds. Then, use the conf_mat() function from dials (part of tidymodels) to create a confusion matrix showing the predicted classes vs. the true classes. Compute the true positive rate (sensitivity), true negative rate (specificity), and accuracy. See this Wikipedia reference if you (like me) tend to forget these definitions. Also keep in mind that a “positive” in this case is a cancellation (those are the 1’s).

preds <- collect_predictions(hotel_lasso_test)
preds

preds %>%
  conf_mat(truth = is_canceled, estimate = .pred_class)

##           Truth
## Prediction     0     1
##          0 34221  7734
##          1  3380 14358

From this confusion matrix:

• The true positive rate will be: \[\frac{34221}{34221+ 7734} = \frac{34221}{41955} = 0.816\]

• The true negative rate will be: \[\frac{14358}{14358 + 3380} = \frac{14358}{17738} = 0.809\]

• The accuracy will be: \[\frac{34221 + 14358}{34221 + 7734 + 3380 + 14358} = 0.8138 \]

• Use the preds dataset you just created to create a density plot of the predicted probabilities of canceling (the variable is called .pred_1), filling by is_canceled. Use an alpha = .5 and color = NA in the geom_density().

preds %>%
  ggplot(aes(x = .pred_1, fill = is_canceled)) + geom_density(alpha = 0.5, color = NA)

Answer these questions:

1. What would this graph look like for a model with an accuracy that was close to 1?

I believe that a model with an accuracy that was close to 1 will have two distinct parts and will not have many overlap points.Moreover, the two parts will be really left-skewed and right-skewed. The reason for that will be: if the dataset said that this person canceled his booking (is_canceled = 1) and the dataset correctly predict that, the prediction will be close to 1. With that, there will have less dataset with .pred_1 < 0.5 while is_canceled = 1 and vice versa for is_canceled = 0. With that, the graph will be split into two parts: one right-skewed red part and one left-skewed blue part.

2. Our predictions are classified as canceled if their predicted probability of canceling is greater than .5. If we wanted to have a high true positive rate, should we make the cutoff for predicted as canceled higher or lower than .5?

As the true positive rate can be calculated as True Positive / (True Positive + False Negative). Therefore, to make the true positive rate higher, we should make the cutoff higher, which then makes the prediction more accurate.

3. What happens to the true negative rate if we try to get a higher true positive rate?

If we try to get a higher true positive rate, we need to increase the cut off rate. As we have less postive cases, we will have more negative cases. With that, we will be more likely to get a false negative cases and the true negative rate will decrease.

8. Let’s say that this model is going to be applied to bookings 14 days in advance of their arrival at each hotel, and someone who works for the hotel will make a phone call to the person who made the booking. During this phone call, they will try to assure that the person will be keeping their reservation or that they will be canceling in which case they can do that now and still have time to fill the room. How should the hotel go about deciding who to call? How could they measure whether it was worth the effort to do the calling? Can you think of another way they might use the model?

For the model, according to the importance variables graph, the hotel should call the person who books the room type P or has been assigned to room I. Additionally, they should call the person who books the room is not refundable.

They could measure who to call based on the rate of that room was being cancelled as well as how likely that will happen. For example, they can look into the variable importance to decide which variables will be important. Once they combine their previous undertanding about the demographic as well as the result from the model, they can find out who should they call.

Another way they can use the model is that they can focus more on some key parts (for example: on the reserved_room P, G or F rather than A,B or C). With these models and its attributions to the cancellation number, they can implement the variables that have positive attributes such as lead_time or stay_in_week_nights for instance to create new rules. For example, they can lower the rate to stay in the hotel during week nights so that they can attract more visitors during that time.

9. How might you go about questioning and evaluating the model in terms of fairness? Are there any questions you would like to ask of the people who collected the data?

In here, since the person who creates the data works in the hotel, one might ask whether he/she collects the data accurately or not. If he/she lowers the cancellation rate during his/her work, it might make the data inaccurate and the model might not be accurate.

Therefore, we should try to see whether to data collecting process might be accurate or not.

Bias and Fairness

Read Chapter 1: The Power Chapter of Data Feminism by Catherine D’Ignazio and Lauren Klein. Write a 4-6 sentence paragraph reflecting on this chapter. As you reflect, you might consider responding to these specific questions. We will also have a discussion about these questions in class on Thursday.

• At the end of the “Matrix of Domination” section, they encourage us to “ask uncomfortable questions: who is doing the work of data science (and who is not)? Whose goals are prioritized in data science (and whose are not)? And who benefits from data science (and who is either overlooked or actively harmed)?” In general, how would you answer these questions? And why are they important?

• Can you think of any examples of missing datasets, like those described in the “Data Science for Whom?” section? Or was there an example there that surprised you?

• How did the examples in the “Data Science with Whose Interests and Goals?” section make you feel? What responsibility do companies have to prevent these things from occurring? Who is to blame?

After reading this chapter, I feel that it is extremely important to ask the question “who is doing the work of data science”. As if the data science jobs consist of looking at the data and trying to find the insights based on that, it will be really important that the data collection part as well as the data analysing process is being done carefully and without any bias. Personally, I have seen some examples of missing datasets, especially when someone tries to manipulate the data for their own purpose.

In my country, sometimes, to prove some certain ideas (such as the literacy rate), they will try to hide or eliminate parts of the data. However, we should ask the question that who is to blame? Is this the person who is inserting the data or is it due to the direction of the upper level class? We will not know the answer to be exact.

All and all, I feel that the best thing that we should do as an individual is try to acknowledge there are still biases in the world today and try the best to prevent it. Moreover, companies need to take steps to make things better and create an environment that will be less biased in the future.